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SUMMARY 


A computer program, designated RETSCP, for the analysis of 
Rocket Engine Thermal Strains with Cyclic Plasticity is 
described in detail. RETSCP is a finite element program 
which employs a three dimensional isoparametric element. 

The program treats elasto-plastic strain cycling including 
the effects of thermal and pressure loads and temperature 
dependent material properties. Theoretical aspects of the 
finite element method are discussed and the program logic 
is described. A RETSCP User’s Manual is presented including 
sample case results. 
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INTRODUCTION 


A new generation of high performance liquid rocket engines 
is being considered for Space Transportation System applic- 
ations. The high performance goal for these engines demands 
high chamber pressures which result in high heat flux levels. 
Engine reusability is a prime objective. With the require- 
ment of thermal and pressure cycling, the stress analyst must 
be able to define the life potential of a given design, 
considering cyclic fatigue where chamber wall stresses are 
sufficiently high to cause plastic strains. 

The state of stress in regeneratively cooled rocket chambers 
varies in three dimensions. For such geometries, a numerical 
method of analysis must be employed. The numerical technique 
which has been given the most attention during the past 
decade is the finite element method. For an outstanding 
introduction to the finite element method, see Zienkiewicz ' s 
text, Reference 1. 

The following report describes a finite element computer 
program designated RETSCP which was developed specifically 
for the purpose of Rocket Engine Thermal Strain analysis 
with Cyclic Plasticity. The program is an outgrowth of a 
General Electric program called ISOPAR, Reference 2. 
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ISOPAR employs a three-dimensional isoparametric element 
to compute the elastic stress distribution in structures 
which can be modeled with relatively few elements. 

The transformation of ISOPAR into RETSCP followed a step- 
by-step approach. First, the program was expanded to allow 
for more elements in the structural model. Then, the 
capability of including thermal loads and computing thermal 
stresses was added. The program was next modified to 
allow non-zero prescribed displacements and to treat sliding 
boundaries. The symmetry condition in a rocket chamber is 
represented by a sliding boundary. Finally, plastic behavior 
with temperature dependent material properties was included. 
In conjunction with this final step, residual strains are 
output on punch cards to allow strain cycle restarts. 

This report begins with a discussion of the theoretical 
aspects of the finite element method. The RETSCP program 
logic and computational scheme are then described. Finally, 
a RETSCP program User's Manual is given which includes sample 
case results. It is intended that a prospective program 
user can go directly to the User's Manual to obtair a 
working knowledge of the program. For application of the 
RETSCP program to specific rocket chamber analyses, see 
Reference 10. 
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FINITE ELEMENT METHOD 


The theory of the finite element method has been well documented 
in several texts (c.f., Reference 1). There are many 
types of elements which have been developel, Reference 3. 

The choice between elements is this: use many simple elements, 

or use few complex elements. The isoparametric element, 
Reference 4, is a very complex element which leads to 
accurate results with a course structural model. 

In this section, the theory of the finite element method 
is described with specific reference to the isoparametric 
element which is used in the RETSCP program. The stress- 
strain analysis, application of boundary conditions, 
thermal loading, and bi-linear plasticity models are discussed 
in the context of the RETSCP program. 



f 


General Theory 

The finite element method is a procedure for approximating 
a continuum by an assembly of distinct elements having a 
finite number of unknowns. For structural analysis, this 
amounts to solving the force-displacement equations for 
the element assembly subject to the prescribed boundary 
values. That is, the folio*. ing system of equations is 
formulated and solved: 

(F) * [K]{6} U) 

where, F and 6 are the forces and displacements at the 
nodal points which connect the elements, and [k] is the master 
stiffness matrix for the assembly. All symbols are defined 
in Appendix A. The appropriate force and displacement 
boundary conditions are used to obtain the solution to 
equation (1). 

The master stiffness matrix is formed by assembling the 
individual stiffness matrices for each element. The element 
stiffness [kj is determined by employing strain energy 
considerations. Apropos to these remarks, the strain 
within each element is related to the element nodal point 
displacements as follows: 

{e> - [B] { 6 } (2) 
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For an elastic structure, the general stress-strain relation- 
ship is 

{o} ■ [D](e) (3) 

Now, the aforementioned energy considerations (c.f, Reference 1) 
imply the following: 

M • /[b] t M [b] dV (4) 

volume 

The functional relationship in equation (2) depends on the 
particular element employed. The detailed manner in which 
the integration, equation (4), is carried out also depends 
on the choice of element. The general procedure, however, 
is to solve the force-displacement equations for the assembly 
under the imposed boundary conditions. 

Isoparametric Element 

Following Reference 1, consider the eight node box element 
shown in Figure 1. The nodal points are located in space 
by their x-y-z coordinates in the rectangular right hand 
system. We introduce a set of parameters (£> n, C) such 
that their values are either +1 or -1 on the element 
faces. A set of eight linear functions of the parameters 
is then defined such that their functional value is +1 
at each corresponding node and zero elsewhere. 
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That is, 

N l » (1/8) (1-5) (1-n) (1-0 

N 2 - (1/8) (1-5) (1 + n) (1-0 

N 3 - (1/8) (1+5) (1+n) (1-0 

N 4 - (1/8) (1+5) (1-n) (1-0 

N £ » d/8) (1-5) (l-n) (1+0 

N 6 - (1/8) (1-5) (1+n) (1+0 

N ? '« (1/8) (1+5) (l+n) d+0 

N g - (1/8) (1+5) (l-n) (1+0 (5) 

Note that these functions apply when the node numbering is 
such that nodes l-2-'-4 go clockwise around the bottom 
when viewed from the top and nodes S-6-7-8 are above 
nodes 1-2-3-4 respectively. 
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Now, the coordinates of any point within the element 

x, y, 2 can be related to the coordinates of the eight 

nodal points x n , y n , z n by the following parametric expressions: 

x ■ N^Xj *• N 2 x 2 + ••* N g x g " {N n } T {x n > 

y - N iyi ♦ N 2 y 2 ♦ ...N 8 y g • {N n } T ly n } 

z ■ N^Zj ♦ N 2 2 2 ♦ ,.,N 8 z 8 ■ {N n } T {2 n } (6) 

Equations (6) thus imply a relationship between (x, y, z) 
and (S,n,; ). 


Bear in mind, that our objective is to evaluate the stiffness 
matrix for the three-dimensional box element, equation (4). 
Thus, we require detailed expressions for the B-m? f 
and D-matrix. The stress matrix, D-matrix, for otropic 

material with elastic modulus E, and Poisson's ra o vis: 


D 


K(l-v) 


(1+v) (l-2v) 


1 v/(l-v) v/(l-v) 0 

v/(l-v) 1 v/(l-v) 0 

v/Cl-v) v/(l-v) 1 0 

0 0 


0 

0 


0 

0 


l-2v 

TIT ^T 
o 
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0 

0 

0 

0 


l-2v 

rrnvr 

0 


0 

0 

0 

0 


l-2v 

TTT^rr 
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The B-matrix relates strai’i at any point in the element 
to the nodal point displacements. The general strain- 
displacement equations are: 
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3w/3x + 3u/3z 
k J 


We relate the displacements of a point in space u, v, w 
to the nodal point displacements (u n ), {v^}, t w n ' as follows: 

u ■ + N 2 u 2 + ...NgUg • { N n > T {u n } 

v - N 1 v l ♦ N 2 v 2 * ...N 8 v 3 - (N„) T {v n } 

w » Njw, ♦ N 2 w 2 + ...NgWg * (N n } T {w n ) (9) 

An element, such as this, for which the same shape function 
expresses the element geometry and displacement fields is 
called an isoparametric element. 
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Substitution of equations (9) into equation (8) gives. 
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To evaluate the displacement derivatives in equation (10), 
we make use of the Jacobian matrix. That is. 
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Substituting equations (b) into equation (11) gives, 
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3N X 

3N 2 . 

..3N S 

TC 
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3N X 
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( 12 ) 


The derivatives in equation (12) are readily obtained by 
differentiating equations (5). This matrix applies for all 
elements and, thus, need only be evaluated once. Then, 
we can determine the Jacobian at any position once the nodal 
point coordinates have been specified. 

It turns out that the derivatives with respect to the 
physical coordinates are related to the parametric 
coordinates as follows: 
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The above matrix defines the elements of the B-matrix 
in equation (10). Thus, upon inverting the Jacobian 
matrix, the B-matrix can be readily evaluated at any 
point in the element. 


Again we restate that our objective is to obtain the stiffness 
matrix, equation (4). Toward this goal we will make use 
of the following relation between element volumes in 
both coordinate systems: 


dV 


xyz “ 


(14) 


where |J| is the determinant of the Jacobian matrix. 

Then, the appropriate form of equation (4) to be evaluated is 
1 1 1 

DO - J J /[bHiOMMI dSdnds (15) 

-1 -1 -1 

Equation (15) is evaluated numerically in the RETSCP 

program. The method employed is wo point Gaussian 

integration based on the following quadrature formula: 

1 

[ f(3c)dx * f (+0 . 57735027) + f (-0. 57735027) (16) 

-1 J 

Of course, the integration is carried out over three 
variables to evaluate equation (15). Thus, the terms in the 
integrand must be evaluated at eight Gauss points within 
the eight node box. 
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One key point remains to be made about the isoparametric 
element used in RETSCP. The element described above was 
based on eight linear shape functions, equations (5). 

The RETSCP element uses those eight functions plus the 
quadratic functions listed below: 

N 9 * 1 -£ 2 

N 1q - 1 -n 2 117) 

Nii * 1 ^ 

Including these, the element has 33 degrees of freedom 
(11 functions times 3 dimensions). Thus, the quadratic 
terms imply a higher oruer element. The functions, 
equations (17), are not associated with any specific point 
in space. For this reason, they are termed nodeless 
variables. The nine internal variables are eliminated 
internally within the program by the technique described 
in Zienkiewicz, Reference 1. Physically this amounts 
to separately minimizing strain energy with respect to the 
variables which are independent of the surroundings 
(otherwise called static condensation, Reference 3). 

Finally, the stiffness matrix is obtained for each isoparametric 
element by the above procedure. Then, the master stiffness 
matrix can be assembled for the entire structure. 
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Boundary Conditions 


Once the master stiffness matrix has been assembled, the 
objective is to solve the governing equations subject 
to the appropriate boundary conditions. That is, to 
solve the system of equations (1), which are rewritten 
below: 
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( 18 ) 


The stress boundary condition is automatically satisfied. 
Namely, forces at nodes on a free-surface are zero in the 
normal direction. 


Prescribed Boundary Forces ; Prescribed force values 
of Pj at the corresponding node are treated simply by 
replacing Fj by Pj in the force vector. 
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Prescribed Displacements : Prescribed displacement conditions 

are treated by modifying the force vector and stiffness 


J 


matrix. Say the jth displacement is to be prescribed 

as . First, leplace F. by F. where 
J J J 

F . * F . - ak . . 

J J J 1 


( 19 ) 


Then, replace the jth row and column in the stiffness matrix 
by zero except kjj which is replaced by 1. This is 
tantamount to eliminating one equation; yet the size 
of the matrix is not reduced. 


As an example of the above procedure, assume Uj has the 
prescribed value a. Then, the resulting equations are 
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( 20 ) 
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Symmetry Condition : The symmetry condition is represented 

by zero displacement normal to the plane of symmetry 
and no restraint along the plane of symmetry (sliding 
boundary) . The symmetry plane is often skew with respect 
to the physical coordinate axis. This is the case for 
a wedge segment with axi-symmetry . Thus, we will derive 
a transformation to treat skew boundary conditions. 

Referring to Figure 2, the displacements in the (x, y) 
system are (u, v) . The skew system (x', y") has a rotation 
of the x-axis of magnitude 6 (positive for rotation of 
x-axis toward y-axis) . The displacements are related 
as follows: 
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The original element properties were evaluated in the 
unprimed system, namely, 

(F> * [K]{6> (22) 
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X 



Figure 2. Notation for coordinate transformation. 


The amount of work done is the same in both systems. 

That is, 

{F'} T {6"} = {F} T { 6 } = {F} T [L]{$'} (23) 

or 

{F' } = [L] T {F} = [L] T [K] [L](6'} (24) 

Thus, we introduce the modified stiffness matrix below 

[r] - [l] t MM (25) 

If, the boundary conditions are introduced in skew coordinate 
directions; then, the corresponding force and displacement 
results are in the skew directions. The entire procedure 
is carried out internally within the program by multiplying 
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the appropriate rows and columns in the master stiffness 
matrix by the appropriate sin-cos terms. It goes without 
saying that only those nodes with skew coordinates need 
be treated. The final results are then transformed back 
into the physical coordinate systems. 

Method of Solution 

The set of governing equations is solved in the RETSCP 
program by Gaussian elimination. The master stiffness 
matrix is partitioned in the interest of computational 
efficiency. The governing equations can be written as 
matrix equations in terms of submatrices. For example, 
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*12 

/ > 
*1 



/ \ 

F ! 

*21 
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F 2 
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\ J 



^ J 


The term Aj is eliminated from equation (26) to give: 

[K*]{A 2 ) = {F*} (27) 

where, 


M * 1*221 ' 

foil 

(28) 

{ F* ) ■ <F 2 ) - 

F 21 ] Fid'hv 

(29) 
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Equation (27) can be solved to give by premultiplying 

by the inverse matrix [K *]" 1 . Then, back substitution 
yields the following: 


u i } * - Ck x1 D Ck 12 ] <a 2 > 


(30) 


Alternately, equation (27) can be partitioned and the 
same procedure reapplied to further reduce the system. 


of including thermal effects is described; also, see 
Reference 5. 
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It should also be noted that 
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lr 

banded matrix. 

This fact also leads to 

a simplification 1 
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in the matrix manipulation. 
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In this section, the method j 
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The temperature difference, referred to a stress free state, 
is input data for each element. Of course, a suitable 
average value must be uced for each entire element. 

The free thermal growth of each element is computed. 

Based on the element stiffness, the nodal forces required 
to mechanically produce the thermal growth are determined. 
These forces are then added to the force vector of the 
entire assembly. Loads and deflections are computed as 
usual for the assembled structure. The stress results 
are adjusted by adding the fully restrained thermal stress 
level for each element. The result is then the actual 
mechanical stress state. 

Bi-Linear Plasticity 

The RETSCP program treats plastic material behavior by 
adjusting the material properties and iterating upon the 
elastic solution. This is the secant modulus procedure 
which was employed in many previous two dimensional finite 
element programs (c.f., References 6 and 7). 

A complete treatment of plastic material behavior is given 
in Reference 8. For the purpose at hand, it is sufficient 
to say that total deformation theory is used; and, yielding 
is based on the Von Mises criteria. For each element in 
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the structure, the average value of the equivalent 
(or effective) stress is computed. That is, the average 
value of the following: 

°e 6 (T 2 x/ T 2 y*’ , 2 xz ) ]' 1 (3i 

Then, according to the Von Mises yield criteria, yielding 

occurs if o is greater than the yield stress from the 
e 

uniaxial stress-strain test. For plastic behavior, equivalent 
stress and plastic strain are related via the uniaxial stress- 
strain curve as shown in Figure 3. 

The RETSCP program employs a bi-linear approximation for 
the uniaxial stress-strain curve. The curve is defined 
by elastic modulus E, yield stress level and plastic 
modulus mE. Plastic modulus and yield can be input as 
functions of temperature. An example of the bi-linear 
stress-strain curve is shown on Figure 4. 


The essence of the secant modulus formulation is as follows. 
First, conduct an elastic structural analysis. Compute 
effective stress and check each element for yielding. 

For elements which indicate yielding, define a new elastic 
modulus called the secant modulus. The secant modulus is 
based on the bi-linear stress-strain curve at the strain 
level corresponding to the elastic result; that is, e tota i. 
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Figure 3 


Relation between equivalent stress and equivalent 
plastic strain. 



Figure 4. Bi-linear stress-strain curve. 




c 


Figure 5. Secant modulus niasticity iteration. 


Associated with e tota i is a bi-linear stress intercept <? new * 
The secant modulus is defined below: 


sec new 
e total 


(33) 


The secant Poisson's ratio, defined to give a consistent 
stress-strain relation, is as follows: 
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sec 


1 1 nE.*,. 

7 ’ly v) sec 


(34) 
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Now, an elastic analysis is again conducted. The stiffness 
matrix, however, is based on E_„ and v co _ for plastic 

ScC 

elements and E and v for elastic elements. The entire 
procedure is repeated and convergence is achieved after a 
few iterative cycles. The process is indicated schematically 
in Figure 5. 

Cyclic Loading 

Two effects of cyclic loading must be considered. First, 
there is the effect of cycling on the material properties 
(see Reference 9). The effect of strain hardening 
(or softening) can be introduced in the program on a 
cycle by cycle basis; or, the cyclic stress-strain curve 
can be input. 

The second effect is the result of plastic deformations 
during one half of the loading cycle. Upon removal of 
the load, residual stresses (or strains) result when 
plastic flow has occured. The residuals, in fact, may 
be sufficiently large to also cause plastic deformation. 

Thus, a stress (or strain cycle) is generated. 


The plastic strain components are related to the stress. 


effective stress, and effective plastic strain as follows; 
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For rocket engine configurations, the shear strains are 
relatively small. Another quantity of interest is the 
equivalent total strain. This value is computed from 
the total strain components as follows: 
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The plastic strain based on this value is then 
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(37) 
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Equivalent total strain, in itself, has no physical significance. 
Within the RETSCP program, the plastic strain components 
and equivalent total strain are computed for each element 
which has yielded. The residual strain components are 
then provided as punch card output for successive run 
calculations. 

The residual strains are read into the program as input 
data for the computation of successive loadings. The 
strains are combined with the thermal strains and analyzed 
in the same manner. That is, the loads at each nodal 
point required to produce the residual strain values are 
computed and added to the assembled load vector. This 
point will be emphasized by example in a later section 
of this document. 
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PROGRAM LOGIC 


The RETSCP program logic is described in this section. 

The general logic is discussed and the program flow iagram 
is given. Some specific points are made concerning the 
subroutine details. The detailed listing of the RI TSCP 
program is given in Appendix B. 

General Logic 

The general RETSCP program logic is to follow the analytical 
procedures outlined in the previous chapter to obtain 
the desired finite element results. 

The computational logic is controlled by the main program 
RETSCP. Subroutines are called as required to perform 
specific calculations. An overlay structure for subroutines 
is employed to reduce core storage requirements. In this 
manner, a specific calculation is performed in a subroutine, 
the results are put onto tape storage (seven tape units 
are utilized), and core storage locations occupied by that 
subroutine are released for reuse. 

The above core storage management procedures allowed the 
RETSCP program size (number of elements) to be greatly 
enlarged from the original ISOPAR program size. In fact, 
the program was enlarged to fully utilize the available 
core storage of the IBM 7094 computer. 


- 28 - 






The data is read into RETSCP from punch cards. For each 
element, the elastic properties and stiffness matrix 
are computed (FEM3D). The master stiffness matrix is formed 
and the boundary values are incorporated (MATRIX). The 
system of equation is solved by Gaussian elimination 
(SOLVE), and the resulting force and displacement values 
at each ncdal ooint are printed out. The elastic stress 
components and equivalent stress values .i" r ' computed for 
each element (STRESS). Now, if the equivalent stress 
exceeds the yield stress a plastic iteration is performed. 

The iteration consists of: first, compute the values of 

secant modulus and Poisson’s ratio (STRESS); then, use 
these values to recompute the elastic properties and stiffness 
matrix for each element (FEM3D) ; finally, complete the 
solution steps above. When the required number of 
iterations have been p rformed, the stress results are 
printed and the residual plastic strains and current 
secant modulus values are punched ou cards to allow cycling 
and restart. 



The flow diagram representing the above steps is given 
in the following section. 
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USER'S MANUAL 


The User’s Manual section contains all instructions necessary 
to prepare data for the RETSCP program. Modeling of the 
structure and preparation of the required input data cards 
are described in detail. Some comments about program output 
are included and sample case results are given. 

Input 

The input for RETSCP consists of punch card data which defines 
the structural geometry, boundary conditions, and materials 
properties. 

The structure is divided into box shaped e’ements which are 
connected by corner nodes. The following procedure for 
locating nodes and elements is quoted directly from Reference 2. 

(a) The 3-dimensional solid is divided by a number 
of non-intersecting surfaces. (Much like slicing a 
loaf of bread.) The surfaces need not be flat or 
parallel, though they frequently are. 

(b) Each such surface is further subdivided by a 
number of non- intersecting lines. (Much like the 
lines on a piece of paper.) The lines need not be 
straight or parallel, though they frequently are. 
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(c) Each such line is further subdivided into a number 
of divisions to give the nodal points. Nodal points 
are numbered in sequence along each line, line by line, 
and surface by surface. 

(d) The nodes on each surface are said to belong to 
the same partition. Partitions are numbered in sequence 
from one side of the solid to the other. (The first 
partition contains the first nodal points.) 

(e) The number of divisions in adjacent lines can 
vary to provide for grading of the mesh. 

(f) 8-noced box elements are formed between adjacent 
surfaces. They are numbered sequentially between each 
pair of adjacent surfaces. The numbering continues 
for successive adjacent surfaces in turn going fmm 

one side to the other of the solid structure. (Although 
in theory the boxes need not be ’’square", it is recommended 
that they be as "square" as the shape of the structure 
permits.) The first element has nodes in the first 
partition. 

The detailed data cards required to execute the RETSCP program 
are listed below. Examples of the data preparation will be 
given in a subsequent section. 
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Card Group 1 ; 


Identification Card 
Number of Cards: 1 

Format: (1114) 


1. Number of partitions (9 maximum) 

2. Number of nodes (225 maximum/25 per partition maximum) 

3. Number of elements (96 maximum/32 per partition maximum) 

4. Number of prescribed displacement nodes (225 maximum) 

5. Number of materials (5 maximum) 

6. Number of degrees of freedom at each node (always 3) 

7. Number of nodes with applied loads (225 maximum) 

8. Starting plasticity iteration number: 1, no iterations 

or 2, iteration starting from elastic solution 
or n, iteration starting from punch card input 
based on iteration number (n-1). 

9. Final plasticity iteration number 

10. Punch output code for successive iterations: 0, no punch 

output 

or 1, provide punch output 

11. Residual stress code: 0, no residual strains input 

or 1, read residual strain card data 
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Card Group 2: 

Coordinate Data 

Number of Cards: 1 per node in order 

Format: 3F16.4 

1. x-coordinate 

2. y-coordinate 

3. z-coordinate 

(inches) 

(inches) 

(inches) 


Card Group 3: 

Node Number Care 
Number of Cards: 1 

Format: 14 

1. Number of nodes 


Card Group 4 : 

Partition Identification 

Number of Cards: 1 per partition in orde* 

Format: 414 

i 

1. First element number in partition 

2. Last element number in partition 

3. First node number in partition 

4. i.ast node number in partition 


( 


• 38 ' 







Card Group 5 : Materials Identification 

Number of Cards: 2 cards per material 

Format: first card 3F16.4 

second card 4F16.4 


Card 1: 

1. 

Young's modulus (psi) 


2. 

Poisson's ratio 


3. 

Coefficient of thermal expansion times 10^(in/in/°F) 

Card 2 : 

1 . 

Yield stress at reference temperature (psi) , Tq 


2. 

Yield temperature gradient (psi/°F) , Aj 


3. 

Plastic modulus times 10^ at reference 
temperature, m Q 


4. 

Plastic modulus temperature gradient 
times 10 6 (1/°F), \ 2 


Note, Card 2 values above are based on the following equations: 


°y = a 0'*l T 


m = mQXlO“3-A2TxlO“6 


(38) 

(39) 


The value of T must correspond to the reference value 
on Card Group 7. 


- 39 - 


% 


M^i 



Card Group 6 : Prestrain Data (can be omitted) 

Number of Cards: 1 per element 

Format: 16* 3F15.8 

1. Element Number 

2. Prestrain in x-direction 

3. Prestrain in y-direction 

4. Prestrain in z-direction 


Card Group 7:' Element Identification 



Number of Cards: 

1 per 

element in order 


Format: 

914, 

F10.3 

1.-8. 

Eight nodal point numbers 



9. 

Material Number 



10. 

Temperature excess over reference 

value 


The eight nodal numbers referred to above are obtained 


for each element: 

(a) Pick a face to be called the top; 

(b) Look down through the top to the bottom face; 

(c) List node numbers clockwise around the bottom 
face (4 values); 

(d) List coincident node numbers clockwise around 
the top face (4 values) starting with the node 
above the first node on the bottom face. 
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Card Group 8 . Element Number Card 

Number of Cards: 1 

Format: 14 


1 . 

Number of elements 



Card Group 9: Displacement Boundary Conditions 


Number of Cards: 1 for each node with 


prescribed displacement 


Format: 414 

, 4F16.8 

1 . 

Nodal number 


2. 

0 if x-displacement is prescribed; 

1 if not 

3. 

0 if y-displacement is prescribed; 

1 if not 

4. 

0 if z-displacement is prescribed; 

1 if not 

5. 

value of x-displacement (inches) 


6. 

value of y-displacement (inches) 


7. 

value of z-displacement (inches) 


8. 

angle of rotation of x-axis toward 

original y-axis (deg.) 


Sufficient displacement boundary condition data must be 
given to fix the body in space. 




r * 

, *T 


i 





Card Group 11 : Iteration Data (can be omitted) 

Number of Cards: 1 per element 


Format : 


16, F20.2, F10.4 


Element number 

2. Secant Young's modulus (psi) 

3. Secant Poisson's ratio 








Output 


The RETSCP output consists of punched cards and printed data. 

Punch cards are provided in conjunction with plastic 
strain analysis. If requested per Card Group 1, the 
secant modulus and secant Poisson’s ratio are punched 
after the final iteration of that run. This allows the 
iterative process to be continued without recomputing 
the initial iterations. For plasticity analysis, the 
residual plastic strain values are automatically punched 
for the final iteration of that run. This data can be 
input directly for subsequent strain cycling calculations 
(Card Group 6) . Secant values and residual strains are 
automatically printed at the end of the printed output 
when the above cards are punched. 

The printed output starts with a list of the input data. 

Note, that the formats may be slightly different from 
the input. For example. Cards 1 and 2 in Group 5 are 
printed in reverse order (Card 2, then Card 1). Also 
Card Groups 3 and 8 are omitted. The input data is 
printed for checking and debug purposes. 


The forces and displacements at each nodal point are 
listed. Values are given in the rotated and rectangular 
coordinate systems. The nodal force data output was 
incorporated to allow numerical evaluation of the net 
section force (such as rocket engine thrust force). 

Detailed stress-strain data is given for each element. 

The stress and strain components at the center of each 
element face are printed as well as the coordinates of the 
face center point. The average stress components for each 
element are also listed. The effective stress which is 
computed in the program is based on the average stress 
components. The yield check data are then summarized 
in the output. This summary consists of effective stress, 
yield stress, total strain, plastic strain, and secant 
values for each element. 

If plasticity iterations are performed; then all of the 
above output data is given for each iteration. Samples 
of output data will be presented as part of the next 
section. 
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Sample Case Results 

Three sample case solutions are presented in this section. 
The cases were selected to demonstrate the capabilities 
of the RETSCP program by successively introducing new 
concepts. Elastic behavior of an isothermal structure 
is treated first. Then, sliding boundaries and plastic 
strains are introduced. Finally, thermal loads and 
strain cycling are illustrated. 

Cantilever Beam : Consider the cantilever beam with concen- 

trated tip load shown in Figure 6. The material is steel 
and the tip load is sufficiently low that elastic behavior 
is guaranteed. The beam is divided into four elements 
as shown in Figure 6. The input data and computer output 
results are presented in Appendix C. The bending stress 
at tne outer fiber is compared with the exact solution 
in Figure 7. The deflection of the nodal points normal 
to the neutral axis (6 Z ) is compared with the exact 
result in Figure 8. This example illustrates that excellent 
results can be achieved with models having few elements. 
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partition number 

Load: P_ * 0.5 lbs. (2.224 Newton) 

Size: L x « 1.0 in. (2.54 cm) 

ly ■ 4.0 in. (10.16 cm) 

L z * 1.0 in. (2.54 cm) 

Matl.: E * 30 x 10 6 psi (20.68 x 10^ N/cm^) 

Figure 6. Cantilever beam sample case configuration 









Thick Wall Cylinder : The second example case is the stress 

distribution in a thick wall cylinder. Due to the symmetry, 
the structure can be modeled by the thin wedge segment 
shown in Figure 9. The boundary condition, with pressure 
load on the inner radius, is zero displacement in the 
tangential direction and freedom to move in the radial 
direction (symmetry condition). The finite element elastic 
stress results for the configuration shown in Figuie 9 
are compared with the exact plane-strain thick wall 
cylinder solution in Figure 10. 

If the stress conditions in the cylinder are sufficiently 

large, yielding will occur. A closed ?orm solution was 

obtained by Mendleson (Reference 8) based on the Tresca 

yield criteria (i.e.,Og -o r >OQ). Yielding under conditions 

of internal pressure will occur from the inner wall to 

some radius p * r ■ p~ . The plastic and elastic stress 

*i 

distributions, according to Reference 8, based on bi-linear 
material behavior are as follows: 

!i ’ h 
°0 P 2 

^0 * C2 
°0 p 2 


Ci(p 2 -1) - P 


Ci(p 2+ 1) + j? 

Or 


♦ C, lnp- p 


+ C jjl + lnp- p | 
' o« / 


p<p c (40) 
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* 2 i 
2 °0 
1 - m 

l-v2m 


f C 2 - m(lV) 

(l-v2m) 

# C 4 - 1 - m 
(l-v2) m 


The quantity 3 is R /Rj and the value of p c is computed from Equation (44) 


below: 
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Stress values for the configuration shown in Figure 9 
were obtained by the finite element method and by closed 
form solution with results shown in Figure 11. 

The difference between the two sets of results is due to 
the different yield criteria employed. Recall that 
RETSCP uses the Von Mises yield criteria; whereas, the 
closed form solution is based on the Tresca criteria. 

Specific input data for the thick wall cylinder example 
is given in Appendix D along with the computed results. 
Note, that the elastic solution is generated as a 
by-product of the plastic analysis (first iteration). 
Summary data only is given for iteration numbers 2, 3, 

4, and 5. 

i 
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Heated Element Cycling ; As a final example we consider 
a single cubic element which is cycled between two temp- 
erature limits. Two opposite faces of the cube are fixed. 
The temperature range is sufficiently great that the 
element stress exceeds the yield stress. Thus, this is 
an example of plastic thermal strain cycling. 

The simple finite element model is shown in Figure 1. 

A sample of the data input and output are given in Appendix 
E. The corresponding stress-strain loop is depicted 
graphically in Figure 13. 

As the material is cooled from its stress free state, 
elastic stresses are built up until the yield point is 
reached (point "a" in Figure 13). Continued cooling causes 
plastic strain to the level indicated by point "b" . The 
total strain at "b" corresponds to the cooling thermal 
strain. The point "c" corresponds to the plastic strain 
residual due to cooling. 

The point "c H is the starting point for the heating cycle. 
As the material is heated, elastic changes occur along the 
line c-d. Plastic changes c<ue to heating occur along the 
line d-e-f. Point "e" corresponds to the residual stress 
state at the original reference temperature. Thus, the 
plastic strain resulting from the cooling half cycle is the 
prescribed displacement for a subsequent analysis. 
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Upon heating the cube, we follow the plastic strain line 
d-e to point "f". The plastic strain at "g M then gives 
rise to the residual stress state "i" as the material 
returns to its original temperature. 

For multi-element structures, the residual stress-strain 
levels during plastic cycling are determined by inputing 
the plastic strain values of all elements and solving 
the residual stress equations for the assembly. 
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Figure 12. Configuration for heated element cycling example. 





Figure 13. Stress-strain loop for heated element cycling 
example. 




APPENDIX A- -SYMBOLS 



F" 

J 

K 

K* 



P 

P 


Matrix of differential functions, Eq. (2), (10) 
Elastic matrix, Eq. (7) 

Modulus of Elasticity 
Secant modulus, Eq. (33) 

Force at nodal point j 
Modified force vector, Eq. (19) 

Equivalent force vector, Eq. (27) 

(in Gaussian elimination method) 

Force vector in skew coordinate system 
Jacobian matrix, Eq. (11) 

Master stiffness matrix, Eq. (1) 

Equivalent stiffness matrix, Eq. (27) 

(in Gaussian elimination method) 

Partitioned stiffness matrix elements 
Element stiffness, Eq. (18) 

Length, or transformation matrix for skew 
coordinate system, Eq. (21) 

Plastic modulus ratio 

Plastic modulus ratio at reference temperature 
times 10^ 

Parametric functions at nodal point n, Eq. (5) 
Load 

Pressure 
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APPENDIX A- -SYMBOLS, Cont'd 


R i 


n 


u 


n 


dV 

v n 

* 

v 


n 
w 

n 

x,y,z 




Inner radius 
Outer radius 
Arbitrary radius 
Radius at yield surface 
Temperature 

Displacement of nodal point n in x-direction 
Displacement of nodal point n in x'-direction 
Differential element of volume 
Displacement of nodal point n in y-direction 

Displacement of nodal point n in y^-direction 

Displacement of nodal point n in z-direction 

Cartesian coordinate system 
Skew coordinate system 
(rotated by angle 0 from x-y) 
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APPENDIX A- -SYMBOLS, Cont'd 


a 


0 

Y xy» Y yz» Y xz 

y P 

T xy 

A 

6 

6" 


"P 

e total 

e et 

e P e P e I 
e x » e y > e z 

e 


v 

v 

o 


sec 


Thermal expansion coefficient 

jth prescribed displacement, Eq. (19) 

Ratio R q /R^, Eq* (44) 

Ratio R /r 
o c 

Shear strains components, Eq. (8) 

Plastic shear strain, Eq. (35) 

Displacement in the partitioned matrix, Eq. (26) 
Displacement matrix, Eq. (1), (2) 

Displacement in the skew coordinate system, Eq. (23) 
Displacement at the nodal point n, Eq. (18) 

Strain matrix, Eq. (2) 

Plastic strain. Fig. 3 
Total strain, Eq. (33) 

Equivalent total strain, Eq. (36) 

Components of plastic 1 strain in x, y, z directions 

Parametric coordinate system. Fig. 1 

Angle of rotation of x-axis into the x'-axis 

in the skew coordinate system 

Poisson's ratio 

Secant Poisson's ratio, Eq. (34) 

Stress 
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APPENDIX A- -SYMBOLS, Cont'd. 


'e 


new 


P 

P, 


l xy 


Effective stress, Eq. (32) 

Yield stress 

New stress, Eq. (33) 

Yield stress at reference temperature, Eq. (40) 
Radial stress component 
Tangential stress component 
Dimensionless ratio r/R. 

l 

r/R A where yield occurs at r 
Shear stress component, Eq. (35) 

Yield temperature gradient 
Plastic modulus temperature gradient 
times 10 6 (1/°F) 


Special Symbols: 


[J. ( ) 

[] T 

ui 

ir 1 


Matrices 

Transposed matrix form 
Determinant value of J matrix 
Inverse matrix 
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APPENDIX B--RETSCP PROGRAM LISTING 
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